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Abstract. Granular media jam into a panoply of metastable states. The way in which these states are 
achieved depends on the nature of local and global constraints on grains; here we investigate this issue by 
means of a non-equilibrium stochastic model of a hindered granular column near its jamming limit. Grains 
feel the constraints of grains above and below them differently, depending on their position. A rich phase 
diagram with four dynamical phases (ballistic, activated, logarithmic and glassy) is revealed. The statistics 
of the jamming time and of the metastable states reached as attractors of the zero-temperature dynamics 
is investigated in each of these phases. Of particular interest is the glassy phase, where intermittency and 
a strong deviation from Edwards' flatness are manifest. 

PACS. 45.70.Vn Granular models of complex systems - 64. 60. My Metastable phases - 45.70.Cc Sandpile 
models - 64.70.Pf Glass transitions 



] 1 Introduction 

Granular media [l] are by now recognised as being pa- 
radigms of complexity [2], especially near their jamming 
' limit [3]. The fact that most grains are too large to be 
. perturbed by the effect of room temperature leads to an 
' 'athermal' dynamics, which is a major cause of this com- 
] plexity - configurations once generated, are remembered, 
■ and their hysteretic effects persist in any ensuing dynam- 
' ics, and the new configurations generated therefrom. An- 
, other origin of complexity in such systems is their generic 
' disorder; this leads in particular to a random landscape 
' of metastable states being explored if a granular system 
I is suitably 'quenched', rather than a crystallisation into 
' an appropriate ordered state. The way in which this land- 
, scape is explored depends strongly on the driving forces 
' applied, in the absence of any real thermodynamics. 

The work we present here is an attempt to explore 
• many of these issues, in particular those to do with the 
nature of ground states in a system near jamming, and 
how these are reached. This is the main motivation for our 
focusing on the effect of 'zero-temperature' dynamics on 
the model we will later introduce, since zero-temperature 
dynamics are known to be the route to systemic ground 
states. Another aspect of interest around which our model 
was designed was the exploration of spatial inhomogenei- 
ties, to answer questions such as: how does position along 
a column of grains influence the dynamics observed there? 
This is an important question for several reasons, one of 
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the earliest being that experimental measurements of den- 
sity along a column of grains revealed wide variations |3j , 
depending on whether they were near the top, the bottom 
or the middle - in its turn, this implied that the com- 
pactivity of the system [51 was non-uniform. A more visible 
manifestation of spatial inhomogcneities is the presence of 
force chains [B] and bridge networks [T\ , which are unique 
signatures of the granular state. 

Given the complexity of issues we wish to investigate, 
we have chosen a minimal model to work with, whose in- 
gredients are based on our experience with several earlier 
ones [819110111] . A feature that all the models share is 
that of orientational disorder; every grain is allowed to 
occupy one of two states, corresponding to 'ordered' and 
'disordered'. Disordered orientations generate voids and 
waste space, whereas ordered ones do not. Implicit in this 
description is the effect of shape, which is most easily un- 
derstood in terms of the rectangular grains of aspect ratio 
a considered in [819] . Grains aligned along their long edges 
(length 1) result in a fully packed column, whereas those 
perched on their short edges (length a) leave voids of size 
I — a. The horizontal orientation is thus ordered, and the 
vertical one disordered. 

Our models became progressively more sophisticated. 
The earliest model of rectangular grains considered in 
[819] was strictly non-interacting, with the only effects in- 
cluded being due to gravity and excluded volume - grains 
could not overlap, and those that were deeper in a col- 
umn were less free to move. One of the major extensions 
in [10|llj was the introduction of the shape parameter e 
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to include grains of arbitrary, i.e., non-rectangular shape - 
such grains can have multiply many orientations in reality, 
but the two-state description was retained in the interests 
of simplicity. The parameter e was allowed to vary over all 
rational and irrational numbers, with a view to describing 
regular and irregular grains, as will be explained further 
below. 

Another difference between the two models was that 
the first one [8 9_ interpolated between jammed and flu- 
idised regimes, whereas the second one was explicitly con- 
structed to examine the jamming regime. In the latter 
case, it seemed reasonable to assume that translational 
diffusion was essentially absent, with all compaction oc- 
curring via orientational rearrangements ~ it was therefore 
appropriate to focus on a column, rather a box [819] . of 
grains. Interactions were introduced in the column model 
of |10|llj such that every grain now not only felt the 
weight of grains above it, but was also constrained by all 
of their orientations. This modelling was appropriate as 
a representation of grains in the free surface layer, which 
were too far from the base to feel an undertow; in the 
jamming limit, however, such grains would be expected to 
feel the orientational constraints arising from grains above 
them. 

The present model incorporates all of the above fea- 
tures, and generalises them to include the presence of ori- 
entational constraints arising from grains below a given 
grain. Unlike our previous models [8|9|10|11] . where in- 
teractions propagated downwards from the top, here they 
propagate both upwards and downwards. Clearly the ex- 
tent of this propagation depends on grain position In this 
sense we model a column of grains with a top, a middle 
and a bottom. 

We devote this paper to a comprehensive investiga- 
tion of the ground states of this model, and how they are 
reached via zero-temperature dynamics. The main ques- 
tions we will answer along the way will be related to sev- 
eral of the issues mentioned above; in particular the issue 
of spatial inhomogeneities will be relevant, since the dy- 
namical regimes attainable via this model will depend on 
which part of the column - top, middle or bottom - is 
being examined. Additionally, we will see that there is a 
panoply of metastable ground states available to the sys- 
tem; the dynamics of their attainment will allow us to 
classify the appropriate regimes as ballistic^ logarithmic, 
activated and glassy. The glassy regime is by far the most 
novel and interesting of these regimes, and will be investi- 
gated in greater depth than the others; its full exploration 
is, however, reserved for future work. 

The plan of this paper is as follows. The definition of 
the model is given in Section 2. Section 3 contains an in- 
vestigation of its static properties, with an emphasis on 
the ground states. Section 4 presents a study of zero- 
temperature dynamics: a rich phase diagram with four 
dynamical phases is revealed and investigated thoroughly. 
A discussion is presented in Section 5. Exact results for 
small systems, as well as other technical results, are pre- 
sented in three appendices. 



2 The model 

Like the unhindered (fully directional) model described in 
[lOlllj , the present model consists of a finite column of N 
grains, labelled by their depth n = 1, . . . , A^. Each grain 
assumes two orientational states, which are referred to as 
ordered and disordered. We set (t„ = -f 1 (resp. cr„ = —1) 
if grain number n is ordered (resp. disordered). A config- 
uration of the column is therefore uniquely defined by the 
orientation variables {(T,i}. There are 2^ such configura- 
tions. 

The model is defined by a stochastic dynamics which 
do not obey detailed balance. In the present context, de- 
tailed balance essentially means a symmetry in m and n 
on the dynamical effect of grain orientation am on (T„ . The 
expressions for the local fields (|2.4|) . (12. 7p clearly do not 
obey such 'action and reaction'. Our model is therefore 
intrinsically out of equilibrium, and its stationary state 
at finite temperature is a genuine non- equilibrium steady 
state 

More precisely, the model is defined as follows. Grains 
are selected in a random sequential fashion and updated 
with the orientation-flipping rates 



w{an = +1 ^ (T„ = -1) = cxp - 



w(cr„ = -1 ^ (T„ = +1) cxp - 



(2.1) 



where, along the lines of previous work [819110111] : 

• -T is a dimensionless vibration intensity, referred to as 
temperature, and related to the 'fast' temperature [T] in 
granular media. 

• An is the activation energy felt by grain n. We make 
the assumption that it increases linearly with the depth 
n, but otherwise does not depend on grain orientations. 
We set 

A„ = — , (2.2) 



Cdyn 



SO that the local frequency 



cxp [ -— ) = cxp 



Cdyn 



(2.3) 



falls off exponentially, with a characteristic length ^dyn- 
This dynamical length corresponds to the depth of the 
boundary layer beyond which grains are frozen out by the 
sheer weight of grains above them. 

• Hn is the local ordering field felt by grain n, which deter- 
mines the orientational response of grain n to the orienta- 
tions {dm} of all the other grains. In previous work [10111] 
the local field Hn only involved the uniform effect of the 
upper grains (m = 1, . . . , n — 1). In the present model, we 
also take into account the back-propagation from grains be- 
low a given grain n (m = n+1, . . . , N). This effect cannot 
be similarly uniform. We assume for simplicity that up- 
ward interactions are exponentially damped, with a char- 
acteristic length ^int, the interaction length. We therefore 
set 

Hn = hn+gjn, (2.4) 
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where we denote by /i„ the uniform effect of grains above 
n (to = l,...,n — 1), and by the non-uniform effect 
of grains below n (to — n + I, . . . , N), whose strength is 
measured by a small (positive) coupling constant g. Both 
components /i„ and j„ of the local field contain the contri- 
bution of every single grain orientation am = ±1 through 
the quantity f{(Jm)- As before |10|llj . the latter assumes 
the following two values: 



/(c^ri 



e 

-1 



if 
if 



-1, 
-1. 



A useful equivalent formula is the following: 



1 



/K) = 2(^-1 -(£ + IK 



In terms of this quantity, /i„ and j„ are: 
ji-i 



(2.5) 



(2.6) 



m— 1 

N 



3n = ^ /(c^m) exp 



m— 71+1 



6. 



(2.7) 



The positive shape parameter e represents an 'effective 
aspect ratio' for a grain of arbitrary shape. This interpre- 
tation originated in the framework of the non-interacting 
model |8l9j with rectangular grains. Rational values of 
e = p/ q imply that the grain size is expressible by a rect- 
angle of sides p and q. Such grains are brick-like, and there- 
fore can be packed perfectly to build some periodic tiling. 
On the other hand, when e is irrational, such tilings can- 
not be built, so that the most close-packed states are those 
of optimal, rather than perfect packing. We thus continue 
to make the following equivalence [10|llj : rational values 
of e imply grains of regular shape, while irrational values 
of e imply grains of irregular shape. 

The parameters of the model are therefore the number 
of grains iV, the coupling constant g, the shape parameter 
e, and most importantly the interaction and dynamical 
lengths ^int and ^dyn- The unhindered model of references 
|9|10|llj is recovered in the absence of coupling [g — 0). 
Throughout the following, we use the notations 



•^dyn 



-l/?dyn 



■^int 



-l/?in 



(2.8) 



To close up, we mention the following recursion rela- 
tions obeyed by the components /i„ and j„: 



hri ~ hn-1 + /((T„-l), 

jn = Xint (/(cTri+l) + jn+l) 



(2.9) 



These relations, together with the boundary values hi = 
jN = 0, provide a fast algorithm to evaluate the local 
fields, to be used in numerical simulations. On the other 
hand, the relations (|2.9[) also imply 



]n-l 



(2.10) 



The latter equation determines the j„ in terms of the /i„: 



''int'tn+l 



] 



(1) 



(1 



_ ,(1) 
Jn I 

N 

E ■ 

i=n+2 



(2.11) 

Throughout the following, we choose to impose for con- 
venience the boundary condition that the uppermost grain 
is ordered: 

(Ti = +1. (2.12) 



3 Statics. Ground states 

The dynamical rules simplify in the zero-temperature limit 
(r ^ 0). Indeed (HH) yields (provided iJ„ ^ 0): 



w{crn = -1 



-1) 



w(a„ 



+1 



-1) 



exp 



00 




if Hn > 0, 
if H„ < 0. 



(3.1) 



Along the lines of references |9|10|llj . a ground state 
of the column is defined as a configuration where the ori- 
entation of every grain is aligned along its local field: 



cr„ = sign Hn = 



if H„ > 0, 
if Hn < 0. 



(3.2) 



The ground states of the unhindered model {g — 0) 
have been investigated in [lOlllj . In that case, the local 
field Hn = hn acting on grain n only depends on the 
grains above n. Equation ()3.2p therefore yields a recursive 
procedure allowing one to construct ground states: 



hn > 
hn < 



= 



+ 1, 
-1, 



hn+1 — hn — 1, 
hn+1 = hn + e. 



(3.3) 



In a ground state, all the local fields hn he in the range 



-l<hn<e. 



(3.4) 



For the present model with a non-zero coupling con- 
stant g, things are more complex. The local field Hn given 
in (|2.4p now depends both on the grains above n (through 
hn) and on the grains below n (through j„)- The condition 
p.2p therefore couples all the orientation degrees of free- 
dom. In particular the ground states admit no recursive 
construction similar to (13. 3p . It is therefore a non-trivial 
task to generate all the ground states of a finite column of 
N grains, whose number depends on g and ^int in general. 

The situation however simplifies in the weak-coupling 
regime [g <C 1), where the ground states can be under- 
stood in terms of those of the unhindered model [g = 0) 
by means of a stability analysis. Just as in the case of the 
unhindered model |10|llj . rational and irrational values 
of e have to be considered separately, in the case of the 
present hindered model. 
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3.1 Rational e (regular grains) 

We first recall some facts about the gromrd state structure 
in the unhindered model. Using (|3.3p in that case, we find 
that whenever the shape parameter e is a rational number 
(in irreducible form) 

e = ^, (3.5) 

q 

the fields /i„ vanish at all depths n such that n — 1 is an 
integer multiple of p + g. The corresponding orientation 
fT„ is left unspecified, and can be chosen at will. Ground 
states are therefore all the random sequences made of two 

well-defined patterns, Pi = H • • and P2 = — + ■ ■ ■ These 

patterns have the same length p + q; they are made of p 
ordered and q disordered grains, and only differ by the ori- 
entations of their two uppermost grains. The model there- 
fore has an exponentially large number of ground states, 
of the form exp{NE), where the configurational entropy 
[12] reads 

r^i^. (3.6) 

p + q 

It turns out that these ground states are still the exact 
ones for the present model, as long as the coupling con- 
stant g is smaller than some threshold gc- In order to show 
this, we approach the problem via a stability analysis. We 
assume that the /i„ field predominates, consider j„ as a 
small perturbation, and see whether the ground states for 
g — are still stable in the presence of the additional 
feedback of the j„ field. 

• Consider first the case when a grain is at a depth n = 
[p -\- q)m + 1 for some integer m = 1, 2, . . . The local field 
hn vanishes, so that the orientation cr„ = 77 = ±1 can be 
chosen at will for g — 0, i.e., in the zeroth order approx- 
imation. Once this choice is made, we have /in+i = /(??), 
a„+i = —77, and hn+2 = e — 1. In order to see if the ground 
state is stable when the effect of the jn field is included, we 
have to first estimate j„ . Equation (|2.1ip can be iterated 
once, yielding 

jn — ^a;int/ln+l + 2;int(l ^ Xint)hn+2 + jn , 

(3.7) 



ji^) ^ (1 _ 2;i„t) 



/ y ^int ^ni + ^int ^ 



m— n+3 



First, we notice that the inequalities 



bound the remainder jn'' as 



''int 



(2) 



allow one to 
< x?.e. Then, 



the values of hn+i and hn+2 obtained above yield the 
inequalities jn > Xint{l — Xint)£ if = -1-1 and jn < 
— a;int(l — 2;int) if <^n = —1- Therefore, the orientation cr„ 
is aligned with the local field Hn = gjn, irrespective of the 
coupling constant g. 

• Consider now a depth n which is not of the form n ~ 
(j3-\-q)rfi~{- 1, so that the local field /i„ is non- vanishing. As 
hn is an integer linear combination of terms /((Tm) equal 
either to —1 or to e p/g, it therefore obeys 



> 



1 



On the other hand, (|2.1ip can again be used to estimate 
The inequalities (|3.4p imply — a;int < jn!^ < a^int^, and 



Jn 



\jn\ < (1 + £)a;int = ^-^Ll x-^^^. 



The inequalities 



and (|3.9p imply 



> 



(j) + g)a;int ' 



(3.9) 



(3.10) 



The full local field i7„ therefore has the same sign as 
hn for all n, provided the coupling constant g is smaller 
than some threshold ^c, so that hn is the dominant field 
in this weak-coupling regime. While the arguments above 
do not allow one to predict exact values of gJB for generic 
rational e, they do yield a lower bound: 



> 



1 



{p + 9)a;int ' 



(3.11) 



Also, similar arguments show that no other ground 
states exist in this model in the same range of allowing 
us to identify, for all g < gc, the ground states of this 
model with those of the unhindered model |10lll| . 



3.2 Irrational e (irregular grains) 

Once again, we review the nature of the ground states 
for irregular grains (irrational e) in the earlier unhindered 
model |10|llj . A unique ground state (corresponding to 
optimal, rather than perfect packing) is obtained for each 
value of e in that case, such that all the fields hn generated 
by the recursion procedure (|3.3p are non-zero. The main 
feature of this ground states is that it is quasiperiodic. 

The nature of the ground states in the hindered model 
under discussion here can be predicted by analogy with 
the low-temperature excitations of the unhindered model. 
Indeed it turns out that the presence of a weak coupling 
{g <C 1) nucleates disorder in this model, in the same 
way as a low but finite temperature (P <C 1) destroys the 
ground state of the unhindered model pT| . 

More precisely, as long as the local field iJ„ has the 
same sign as /i„, i.e., for |/i„| > g\jn\, the ground states of 
both unhindered and hindered models are identical. The 
depth up to which this stability condition is satisfied can 
be estimated as follows. Equation (|3.9p shows that typi- 
cal values of jn are of order ccint . The stability of a given 
ground state is determined by the grains n where hn and 
jn are comparable, i.e., |/i„| ~ gx-mt ^ 1- Such sites with 
very small hn fields are nothing but the nucleation sites 
which dominate the low-temperature behaviour of the un- 
hindered model [11]. The typical distance L{g) between 
two consecutive nucleation sites diverges as 



Ha) 



1 



gXint 



(3.12) 



^ The exact threshold coupling gc will, however, be deter- 
mined later on for e = 1 (see ^^). 
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in the regime of a weak coupling (gxi^t *C 1). Thus - as 
one might expect - the closer g is to zero, the less chance 
there is that disorder is nucleated. 

As long as the size TV of the column is smaller than 
h{g), or, equivalently, the coupling constant g is smaller 
than l/(A^a:int), the unique ground state of the model is 
the quasiperiodic ground state of the unhindered model. 
For larger columns, a typical ground state can be thought 
of as a sequence of independent quasiperiodic patches of 
mean length h{g), pasted together end on end. For very 
large sizes (TV ^ ^(g)), there is therefore an exponentially 
large number of such ground states. The corresponding 
configurational entropy: 



1 



L(5) 



gXint 



(3.13) 



vanishes in the weak-coupling regime {gxint ^1). 

Finally, we can provide an integrated description of the 
ground states of an irrational s and of its rational approx- 
imants. For g small and e a fixed irrational, we consider 
the sequence of its rational approximants Sk = Pk/lk [H^- 
The periods pk + Qk of these approximants typically grow 
exponentially fast with the approximant order k. Equa- 
tion (|3.1ip implies the following. For the first rational ap- 
proximants in the series, whose periods are smaller than 
l/{gxint), the ground states are the same as in the un- 
hindered model. For all the higher approximants, whose 
periods exceed this threshold, ground states are expected 
to be made of nearly independent patches, whose char- 
acteristic length L{g) is given by (|3.12[) . just as for the 
limiting irrational. Notice that the period of the approx- 
imant that divides these two behaviours is of the order 
of L((7) i/{gxint) - there is a reassuring consistency in 
this. 



3.3 More details for e = 1 

We now focus on the simplest case, which is obtained when 
the shape parameter equals s = 1. The complex behaviour 
we obtain even from this simplest of all cases is a tribute 
to the inherent richness of the model. 



Generic configurations 

Assume that the column size N is even for definiteness. 
Consider first a generic configuration. Equation (|2.5p sim- 
phfies for e = 1 to /(cr„) = — cr„. As a consequence, the 
components of the local fields read 



m— 1 



N 

E ^ 

m— n+l 



(3.14) 



The first expression shows that hn is an integer, whose 
parity is opposite to n. Therefore: 

• When the depth n = 2fc is even, hn is odd, and thus 
always non- vanishing. In the weak-coupling regime {g <C 
1), we thus obtain iJ„ w /i„, irrespective of ^mt- 



• When the depth n = 2A: — 1 is odd, /i„ is even, and may 
therefore vanish. When it does, we have Hn = gjn, so that 
the sign of Hn depends on ^int in general. 

In order to understand the complexity which can arise 
from this dependence on the interaction length ^int, let 
us focus on a column of size iV = 6, in the particular 

configuration H 1 h- f. We have as = -1-1, /13 = 0, 

and H3 = gj3 = ga;int(l — a;int — xf^^). The parenthesis 
is a second-degree polynomial in x-mt- It vanishes (in the 
physical range < Xjnt < 1) when x-mt equals the inverse 
golden mean 



%/5 - 1 



0.61803. 



(3.15) 



Thus when < x-^t < 4>, H^ = gj^ is positive, and the 
condition CT3 — signiJa is fulfilled. This condition does not 
hold for <j) < Xint < 1- 

This example is illustrative of a general property of the 
model. For a finite column made of N grains, the statics 
and dynamics of the model depend on the relative position 
of Xint with respect to a finite number of threshold values, 
where quantities of interest are discontinuous in general. 
These threshold values are given by the following rule: 
they are all the roots (in the physical range < Xi^t < 1) 
of all the reduced fields jn/xint with n > 3 odd, viewed as 
polynomials in a^int, in all the configurations. The polyno- 
mial jn(a^int)/a^int has even degree d — N — n — 1, so that 
d<N-4:. 

These threshold values all lie in the range 1/2 < Xint < 
1. This fact can be seen as follows. The expression p.l4p 
for the local field jn can be recast as 



Jr, 



-XintfJ; 



N 

X\y^+ Or] 



n+l / , -^int 

m— n+2 



(3.16) 



The sum in the right-hand side is smaller than x?j^/(l — 
Xint) in absolute value. As a consequence, the local field j„ 
always has the sign of the leading term, and therefore can- 
not vanish, as long as Xint > 2^hit/ {^ — Xint)i i-e., x-^at < 1/2. 
The above property is responsible for a strikingly general 
result: no dynamical quantities for arbitrary system sizes 
depend on xi^t for < Xint < 1/2, i.e., for ^int < 6nt,i, 
with 

6nt.i = 7^~ 1.44269. (3.17) 
m 2 

In particular, this explains the plateau in the data to be 
presented in Figure [31 

We have generated all the relevant polynomials up to 
degree d = 12. The numbers Ad of physical roots with 
degree d, and the numbers Bat of threshold values for a 
column of N grains, are found to be the foUowingQ 



A2 = 1, A4 = 5, Ae = 19, As = 97, 
Aw = 442, A12 = 1880, 
-B2 = 0, B4 = 0, Bq = 1, Bs — 6, Bio 
B12 = 121, Bu = 563, S16 = 2443. 



25, 



(3.18) 



^ Since some of the polynomials are reducible, the same root 
can be generated several times, so that the number Bn of 
distinct thresholds may be smaller than the sum A2 + A4 + 

\-An-4. 
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We can extract several insights from the above numbers. 
Since B2 — B4 = 0, columns of size N = 2 and = 4 
exhibit no dependence on ^int at all. For = 6, we 
have Bq = I: the unique threshold value is the inverse 
golden mean introduced in p.l5p . Consequently, quan- 
tities typically assume two different forms in the intervals 
< Xint < (j) 01 (j) < Xint < 1- For = 8, we have Bs = 6: 
quantities typically assume seven different forms in the in- 
tervals demarcated by the six threshold values of xint , and 
so on. The interested reader is referred to Appendix A for 
the explicit verification of these predictions for columns of 
size A^ = 4 and A^ = 6. 



Ground states 

We now focus on the ground states in the simple case 
where e = 1. They consist of dimerised configurations, 
made of the patterns Pi — ^ — and P2 = — h . Assuming 
again that the column depth A^ is even for definiteness, 
the generic ground state can be described as follows: 



with 



0'2fe-l = 



= -<J2k=Vk {k = l,...,N/2), 



(3.19) 



with the dimer variable rjk 
responding to the Pi (resp. 
condition (|2.12p implies r/i : 
ground state read as follows: 



= +1 (resp. rjk = — 1) cor- 
P2) pattern. The boundary 
= -1-1. The local fields in a 



N 



N 



£1 = — ^ jnCTn- 



(3.24) 



This definition can be motivated as follows. If the cr„ were 
independent spins in external fields H„, (|3.22|) would be 
the corresponding Hamiltonian. In the present model, we 
recall that the local fields i/„ depend on the orientations 
<7rn in a complex and non-symmetric way, so that the dy- 
namics does not obey detailed balance, and the statics 
is not described by any simple Hamiltonian. The pseudo- 
energy defined by (j3.22|) . however, provides a useful mea- 
sure of the amount of configurational disorder in the full 
interacting system and, in particular, allows us to classify 
the ground states. 

In the case of a ground state, the first component of 
the pseudo-energy reads Eq ~ —N/2 (irrespective of the 
ground state), whereas the second one reads 



£1 = - 



NXint 



(1-Xi„t)' 



„2l~2k-l 



l<k<l<N 



VkVh (3.25) 



in terms of the dimer variables 77^. The first term in the 
above expression, 



Nx;u 



(3.26) 



^2fe-l = 0, h2k = -Vk, 



N/2 



J2fe-1 = XintVk - (1 - Xint) ^ 

l=k+l 

N/2 

32k = -(1-Xint) 2^ Xj^t m- 

l=k+l 



2l-2k 



(3.20) 



There are 2^/^ ground states in total, or 2^/2^1 if ((TT^ 
is taken into account. 

The two crystalline (uniform) ground states 



(3.21) 



play a special role; intuitively, this is because any other 
ground state has conflicts between successive dimer pairs, 

for example in a configuration such as H 1- where the 

second and third orientations would be in non-ideal posi- 
tions. Such conflicts would of course also be present, albeit 
more weakly, between dimer pairs that were more distant 
from each other along the column, e.g. -I — • • • — h. 

This observation can be turned to a quantitative clas- 
siflcation of ground states by means of the pseudo-energy. 
This quantity is defined for an arbitrary configuration as 
follows: 

N 

£ = -Vff„a„, (3.22) 



I.e., 



f = ^0 + .9^1, 



is the mean of £1, in the sense of a uniform average over 
all the ground states. The second term in (|3.25p represents 
the fiuctuation in £1 from a ground state to another, which 
typically grows as N^^'^. 

The crystalline states U± introduced in p.2ip . respec- 
tively corresponding to rjk = +1 and 77^ — —1 for all k, 
are the two absolute (global) minima of the pseudo-energy. 
We find 



Nxu 



a;int(l 



'^int) 



l + a:int (H-a;int) 



(3.27) 



It follows that the crystalline ground states are separated 
from a bulk of roughly equivalent metastable states by an 
extensive pseudo-energy gap 



A£ ^ {£1) - £i{U±) 



NXjntjl - Xint) 
2(1 -I- Xint) 



(3.28) 



It is remarkable that our model generates the kind of (free) 
energy landscape that is familiar in theories of glasses [13] , 
and guessed to be valid for grains in the jamming limit [1], 
where crystalline states lie well below a band of metastable 
states. 

Finally, the exact threshold coupling gc for ground- 
state stability for e — 1 can be evaluated as follows. Equa- 
tion (|3.20p implies that the state U+ is the first to be 
destabilised by an increase of the coupling, and that its 
weakest link is the second grain orientation (72 . The corre- 
sponding local field reads H2 = —l+gj2, where j2 assumes 
(3.23) its largest possible value j^"""^ = Xint{l - x^^'^j/{l + Xint)- 



J.M. Luck, Anita Mehta: Dynamical diversity and metastability in a hindered granular column near jamming 



7 



The threshold coupUng constant at which this first desta- 
bihsation takes place therefore reads gc = the 
limit of an infinite column, we therefore obtain the follow- 
ing exact expression for the threshold coupling: 

= = ei/«'- + 1. (3.29) 

•^int 

The threshold coupling is found to be a decreasing func- 
tion of the interaction length ^int, blowing up exponen- 
tially at small ^int, and reaching a finite asymptotic value 
2 in the ^int oo limit. It is interesting to observe that 
the general bound (|3.1ip . i.e., gc > e^/'^'"'/2 in the present 
case (p — q = 1), captures the qualitative features of the 
dependence of the exact result on ^i^t. 



4 Zero-temperature dynamics 

The application of zero-temperature dynamics is the can- 
onical way of finding the ground states of a model. In the 
present model, grains are aligned one at a time with their 
local fields. More precisely, grain n is selected at a rate 
given by (|2.3|) . and its orientation variable cr„ is aligned 
along the local field i/„ introduced in (|2.4p . according to 
the deterministic rule 

cr„^ sign ff„. (4.1) 

This rule is well-defined for a non-zero coupling constant 
g, because the local fields i7„ generically do not vanish. 

This leads to metastability in the following sense: a 
finite column of N grains in an arbitrary initial config- 
uration is eventually driven to an absorbing configura- 
tion or attractor, in a finite jamming time T. This at- 
tractor is necessarily one of the ground states described 
earlier, i.e., a configuration where every orientation cr„ 
is aligned with i/„. For a coupling constant g less than 
the threshold gc given in (|3.29p . the attractors are the 
2^/2 dimerised configurations, whose number is halved to 
2N/2-1 j£ ^YiQ boundary condition (I2.12p is taken into ac- 
count. Arbitrary initial conditions can lead to any one of 
these metastable configurations being reached; they are 
however fragile in the sense that a slightly different ini- 
tial condition or stochastic history generically leads to 
another attractor being reached instead. This fragility of 
metastable states is one of the characteristics of granular 
media [1]. 

In what follows, we will focus on two aspects of zero- 
temperature dynamics: 

• Statistics of the jamming time. The jamming time T is 
the random time the system takes to converge to an at- 
tractor, it being understood that the initial configuration 
is disordered and randomly chosen. The N dependence of 
both the mean jamming time (T) and of its full probability 
distribution are of interest. In this respect, we introduce 
for further reference the reduced variance 



• Statistics of the attractors. The statistics of the attrac- 
tors reached by stochastic dynamics is also of special inter- 
est, especially in relation to Edwards' flatness hypothesis 
[5]. We anticipate non-trivial results only in the glassy 
phase, which will be investigated in Section 4.4. 

Inspired by previous work |10lll| , we monitor the vari- 
ous dynamical regimes of the model by means of the thick- 
ness L(t) of the upper ordered layer of the column, defined 
as the depth of the uppermost grain which is not aligned 
with its local field: 

L(t) = inf{n|cr„^signiJ„}. (4.3) 

Figure [1] shows typical tracks L{t) for four representa- 
tive choices of values of ^int and ^dyn- Here and throughout 
the following, we choose for definiteness the value g = 0.01 
of the coupling constant. This value is deep in the weak- 
coupling regime, where all the results are virtually inde- 
pendent of the precise choice made for the coupling con- 
stant The dynamical behaviour observed turns out to 
be very strongly dependent on the lengths ^int and ^dyn, 
and suggests the existence of four qualitatively different 
dynamical phases, which we have named ballistic, loga- 
rithmic, activated and glassy. They will be investigated in 
greater detail in what follows. The dynamical phase dia- 
gram in the ^int — ■Cdyn plane presented in Figure [5] shows 
already the existence of genuine phase boundaries (where 
crossover phenomena become arbitrarily sharp in the limit 
of an infinite column) denoted by full lines, and crossover 
phenomena (which occur when ^dyn is comparable with 
the column size N) indicated by dashed lines. 

4.1 Phase I (Ballistic) 

This phase is illustrated by panel I of Figure [1]- it corre- 
sponds to small ^int and large ^dyn- Physically this implies 
that one is looking at layers near the free surface (large 
Cdyn) of a column where grains feel correlations from below 
only weakly (small ^mt)- In other words, we are looking at 
the 'top' of a granular column. 

The thickness L{t) is observed to grow on average lin- 
early with time: 

im) « Vt. (4.4) 

The (relatively small) fluctuations between different tracks 
correspond to different stochastic histories. Equation (|4.4p 
shows that an ordered layer propagates ballistically down 
the column with velocity V, a phenomenon which was al- 
ready encountered in the unhindered model [10|llj . When 
L{t) becomes equal to the column depth N, an attractor 
is reached and the dynamics stops, so that 



again up to relatively negligible fluctuations. 

Figure [3] shows numerical values of the velocity V, ob- 
tained by averaging L{t) over many independent initial 

For comparison, we recall that the threshold coupling Qc 
(see p.29p l is always larger than 2. 
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100 200 300 

t (10^) 

Fig. 1. Plots of the thickness L{t) of the upper ordered layer 
of the column against time t, for g — 0.01, illustrating the 
four phases. I: Ballistic phase (several tracks with ^int = 3, 
Cdyn = oo, = 500). The slope V = 6.95 of the thick line is 
taken from Figure|3] II: Activated phase (one single track with 
^int = 50, 4dyn ~ oo, TV — 70). IIL Logarithmic phase (several 
tracks with ^int = 3, Cdyn = 200, A = 500). The thick line 
shows the result (|4.27|l with V = 6.95. IV: Glassy phase (one 
single track with ^int = 50, .fdyn — 7, N = 70). 



configurations and histories (at least 10^ per point). The 
first interesting feature is a plateau region observed for 



00 
A N 

t 



I (Ballistic) 


II (Activated) 




III (Logarithmic) 


I'V^jciassy)^^ 







1/a 



00 



int ' 

Fig. 2. Schematic zero-temperature dynamical phase diagram 
of the model in the Cint— Cdyn plane, showing the four dynamical 
phases revealed and investigated below. The numbers ^int,c and 
floo will be determined later (see (|4.1ip and (|4.25|l ). 
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Fig. 3. Plot of the ballistic velocity V against ^int for ^dyn ~ 
oo. Arrows indicate the value ^int,i (see (|3.17p ) below which 
V is constant and equal to Vo given in (|4.6p . and the critical 
point ^int,c (see (|4.1ip ') at which V vanishes linearly. 



^int smaller than the threshold value ^int,i determined in 
p.lTp . As predicted in Section 3.3, the velocity V is found 
to be constant over this region, and equal to its value in 
the ^int limit: 

Vo « 9.75. (4.6) 

As ^int increases, the effects of frustration begin to kick 
in more and more via the back-reaction j„; the resulting 
inefficiency of the zero-temperature dynamics causes the 
velocity to decrease progressively with fint- 

It seems strange that Vq is such a large number, espe- 
cially given that it leads to the apparition of other large 
dimensionless numbers, such as ^int,c ~ 28.4 (see (|4.1ip ) 
or Dc ~ 37 (see (I4.15|) ). Fortunately, we can provide a 
simple explanation for the high value of Vq , that the nat- 
uralness principle IJA^ would demand. We recall (see the 
lines following (|3.16p ) that the component j„ of the lo- 
cal field always takes the sign of — ct„+i for ^int < Cint,i- 
Assume that the uppermost 2k grains of the column are 
already dimerised. So far as zero-temperature dynamics 
is concerned, the next two orientations a2k+i and a2k+2, 
are entirely decoupled from the rest of the column. Indeed, 
h2k+i = 0, so that H2k+i = 9j2k+i has the sign of -a2k+2, 
whereas H2k+2 ~ ^2fc+2 — ~cr2fc+i- If we make the sim- 
ple assumption that the four configurations of these two 
orientations are equally probable, the mean time it takes 
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to go to one of the two possible attractors can be shown 
to be 1/4. This newly formed dimer has a spatial size 2, 
and the corresponding front velocity is 2/(1/4) = 8. The 
actual velocity Vq of (|4^p is only about 21% above this 
simple-minded estimateQ. 

The data of Figure [3] also show that V vanishes linearly 
as the borderline between Phases I and II is approached, 
i.e., as fint — > Cint c This behaviour fits in a natural way 
within the effective description put forward in Section 4.2 
in terms of biased Brownian motion (see (|4.19p ). 



4.2 Phase II (Activated) 

This phase of relatively large ^int and ^dyn is illustrated by 
panel II of FigureHJ One is still considering grains that are 
relatively free to move, as in the top layers of a column, 
but now grains are increasingly constrained as a result of 
grain orientations below them. 

There is, typically, only very weak order in the col- 
umn before it happens to jam: this is exemplified by a 
mean thickness {L{t)) which is quite small compared to 
the column depth N. Also, L{t) exhibits wild fluctuations 
around its mean, which look stationary over the very long 
time it takes for the system to jam. After sporadic excur- 
sions to larger values, the layer thickness suddenly jumps 
to L{t) — N, so that an attractor is reached. 

This phenomenology is typical of an activated phe- 
nomenon. We therefore expect that: 

• The statistics of the jamming time T should be approx- 
imately given by an exponential distribution: 



p(r) = — exp — 



(4.7) 



characterised by a single scale (T) , with unit reduced vari- 
ance {Kt = 1)- 

• The mean jamming time should grow exponentially with 
the column size: 



(4.8) 



at least for very large N, where a(^int) is the reduced ac- 
tivation energy per grain, i.e., the height of the entropic 
barrier the system has to cross to reach the ground state. 

Huge finite-size effects rule out an accurate numeri- 
cal exploration of the activated phase for generic ^int and 
^dyn- lu the following, setting ^dyn ~* oo, we examine two 
limits of particular interest. We explore first the crossover 
between Phases I (ballistic) and II (activated), as ^i„t ap- 
proaches the critical value ^int.c- Next, we examine the 
regime of deep activation, when the effect of upward in- 
teractions is maximal (^int ^ N). 



^ We mention for comparison that in the analogous case in 
the model of [10111] . that of an irrational e — > 1* in the im- 
mediate neighbourhood of £ = 1, a similar line of reasoning 
yields the value 2 for the velocity, whereas the measured ve- 
locity Vji « 2.38 is about 19% above that estimate. 
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Fig. 4. Top: logarithmic plot of the ratio {T)/N against ^int, 
for a dynamical length ^dyn = oo, and variable column size A*'. 
Bottom: plot of the reduced variance Kt against ^int, for the 
same parameters. 



Crossover between Phases I and II 

In order to understand the crossover between Phases I and 
II, we invoke the following picture of the behaviour of the 
thickness L{t) of the ordered layer. In either of the two 
phases, it starts from the surface and eventually propa- 
gates to the base, when an attractor is reached. In the 
purely ballistic case (^int very small) , the layer essentially 
shoots down to form an attractor. The effect of increas- 
ing 6nt is to 'admit impediments' to this pure flow, to 
cause Lit) to fluctuate (diffuse) increasingly before the 
whole column reaches an attractor. The uniform effect of 
grains above any given grain and the back-reaction of the 
grains below it are responsible for the frustration that is 
increasingly encountered in the search for an attractor as 
^int increases. The value of ^int at which both interactions 
balance out is the critical point ^int,c, at which the velocity 
V vanishes, so that the dynamics is purely diffusive. 

The above intuitive picture is corroborated by Fig- 
ure 21 showing a logarithmic plot of the ratio {T)/N and 
a plot of the reduced variance Kt against ^int, for sev- 
eral values of the column size N. There is evidence of a 
continuous phase transition between a ballistic phase for 
6nt < 6nt,c and an activated phase for ^int > Cint,c- We 
notice from the top panel that (T) /N is roughly indepen- 
dent of N in the ballistic phase, whereas it grows fast with 
N in the activated phase. On the other hand, the plots of 
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Fig. 5. Top: logarithmic plot of {T)/N'^ against iV(Cint-Ciiit,c). 
for a dynamical length ^dyn = oo, ^int.c = 28.4, and variable 
TV. Bottom: plot of the reduced variance Kt against X, for 
the same parameters. Full lines: plots of the analytical results 
(|B.31|) and (|B.32|) for the effective model, rescaled according 
to (f4l2ll . 



(jB.32p . derived analytically in Appendix B for the effec- 
tive model. The excellent quality of the accord suggests 
that the effective model predicts the exact finite-size scal- 
ing functions of the model. The best agreement, shown 
as full lines in Figure [51 is obtained with the following 
identification: 



X 



-115(z + 3.4), ln(To/iV^) 



-4.3 



(4.12) 



between, on the empirical side, the scaling variable X of 
the column model introduced in ()4.10|) . and, on the theo- 
retical side, the scaling variable 



VN 

15' 



(4.13) 



of the effective model, introduced in (|B.30|) . and the dif- 
fusive time scale 



2D' 



(4.14) 



introduced in (|B.24p . This last relation, together with the 
second equation of (|4.12p , enable us to predict the critical 
value Dc of the diffusion constant for ^int ~ 6nt,c- We thus 



obtain Dr 



,4.3 



/2, i.e. 



D, » 37. 



(4.15) 



The reduced variance of the jamming time at the critical 
point is predicted in (jB.26|) to be a universal number: 



(4.16) 



Kt approximately cross at a critical value Kt ^ 0.7 for 

Cint = Cint,c 26. 

This picture can be turned into the following effective 
model. We treat the thickness L{t) as a collective coordi- 
nate, and model its dynamics by a biased Brownian mo- 
tion on an interval, with velocity V and diffusion constant 
D. The motion starts at time i = at an initial point 
L(0) very near the free surface of the column, which is 
considered as a reflecting boundary. It ends at the ran- 
dom hitting time t = T when L{t) visits the base of the 
column, i.e., L{t) — N, for the first time. Accordingly, the 
base is considered as an absorbing boundary. This effective 
model is analysed in detail in Appendix B. 

Figure [5] shows that both the mean jamming time (T) 
and its reduced variance Kt obey finite-size scaling laws 
of the form 



{T)kN^F{X), Kt^G{X), 



with 



^ = A^(6nt-?i„t,c). 

The best data collapse is obtained for 
6„t,c ~ 28.4. 



(4.9) 
(4.10) 

(4.11) 



Furthermore, the data are observed to be in accurate 
agreement with the finite-size scaling results (|B.31|) and 



This prediction is again in good agreement with the ap- 
parent crossing point of the data of the lower panel of 
Figure HI 

Finally, we compare the predictions of the effective 
model in the ballistic and activated phases with the above 
results. 

• Toward the ballistic phase (^-mt < Cint.c, *-e-, V > 0). 
In the ballistic phase, the prediction (jB.2ip for the mean 
jamming time: 

N D 

iT)^y^^ (4.17) 

exhibits the observed ballistic behaviour (|4.5p . up to a 
finite negative correction due to diffusion. 

The fluctuations of the jamming time around its mean 
are predicted to be Gaussian, with a reduced variance 
given by (|R22| : 

Kt w — . (4.18) 

This fall-off as 1/iV agrees with the observation made 
above that fluctuations become relatively negligible for 
large columns. 

The critical regime of the ballistic phase corresponds 
to the regime X — > — oo, i.e., z — s- +oo, in the finite-size 
scaling laws (|4J|) . Equations (IB.3ip and (|4.12|) imply that 
the scaling function F{X) falls off as F{X) « Af/\X\, 
with Ap w 115/D w 3.1. This estimate implies in turn 
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that the velocity vanishes linearly as ^int 



^>int,c' 



V « 0.32te, 



6nt)- 



as V : 



(4.19) 



The numerical value of the prefactor is in good agree- 
ment with the slope of the extrapolation curve shown as 
a dashed line in Figure [31 

• Toward the activated phase fi^int > Cint,c; V < 0). 
In the activated phase, the prediction (jB.27|l for the mean 
jamming time: 

D 



(T) 



\V\N/D 



grows exponentially with N ^ as anticipated in (j4 
corresponding activation energy per unit length 

1^1 



(4.20) 
The 

(4.21) 



is essentially given by the negative of the velocity. The 
reduced variance of the jamming time predicted by (|B.29p : 



K^ 



1 



2{VN + iD) ^^\v\N/D 



D 



(4.22) 



converges exponentially fast to its limiting value unity, 
characteristic of an exponential distribution. 

The critical regime of the activated phase corresponds 
to the regime X +oo, i.e., z — > —oo, in the finite- 
scaling laws ()4.9p . The effective model predicts that the 
scaling function G has an exponential convergence toward 
G'(-l-oo) = 1, whereas the scaling function F grows ex- 
ponentially as F ~ exp(— z) ~ eyi-p{BpX) with Bp ~ 
1/115 ~ 0.0087. These results corroborate our expecta- 
tions, including (|4.7[) and (|4.8[) . They also imply that the 
activation energy per grain vanishes linearly as ^int ~^ 
?int c as a(^i„t) « Bp{S,,^t - Cint,c), i.e.. 



a(6nt)« 0.0087(Cint-6nt,c). 



(4.23) 



Before leaving this topic, we emphasise that the simple 
picture of a Brownian particle whose velocity V changes 
sign at ^int,c appears to explain all our observations on 
this crossover. 



Limiting behaviour for ^int ^ N 

We now look at the slowest possible dynamics in the acti- 
vated phase; this will clearly occur when the column is at 
its most correlated, where ^i„t is much larger than N , still 
keeping ^dyn = oo for simplicity. We are thus led to inves- 
tigate the doubly singular limit where ^dyn = 6nt = oo- 
The only free parameter is then the column size A'^. 

Figure [6] shows plots of the mean jamming time (T) 
and of its reduced variance Kt against N . When the col- 
umn size N is rather small, the system still behaves more 
or less ballistically; this fast dynamics leads to the nearly 
linear growth of (T) with iV, and the concomitant de- 
crease of the variance Kt as a function of N . When N 
is large enough, the system is fully activated, so that the 
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Fig. 6. Top: logarithmic plot of the mean jamming time (T) 
against A'', for ^int = Cdyn = oo. Full line (hardly visible): fit 
(T) = (In 2)/2 - 6.515 In A + 15.04 to the data for A > ^int.c 
Bottom: plot of the reduced variance Kt against A, for the 
same parameters. Arrows show the crossover scale A = ^int.c 
(see dtni). 



jamming time grows exponentially with N , while the vari- 
ance increases, rapidly converging to its asymptotic value 
Kt = 1. The crossover between these behaviours occurs 
when the column size N is of the order of ^int.c (shown as 
arrows in both plots). 

Finally, we focus on the activation energy a(^int) de- 
fined in (|4.8p . When ^int — > oo, the column is at its most 
activated; hitting an attractor is then a totally random 
process. The jamming time is therefore expected to be sim- 
ply given by the ratio between the number Qq = 2^ of dis- 
ordered initial configurations and the number flea = 2^^"^ 
of possible attractors: 



(T> 



2N/2 



(4.24) 



A similar purely entropic result is shown in Appendix C 
to hold within a toy model of a Markovian dynamics on 
an assembly of independent two-level systems. The result 
(j4.24p implies that a(^int) saturates to the value 



In 2 
~2~ 



0.34657. 



(4.25) 



This limiting value of the activation energy has been in- 
corporated into the fit presented in the upper panel of 
Figure [51 The good quality of the fit can be viewed as 
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Fig. 7. Plot of (1/iV) ln(r) against l/^dyn for iV = 50. Empty 
symbols (^int = 3) demonstrate the crossover between Phases I 
(ballistic) and III (logarithmic). Dashed line: prediction (|4.28p . 
with V — 6.95. Full symbols (^int ~ oo) corroborate the 
crossover (|4.33|) between Phases II (activated) and IV (glassy) , 
in spite of large finite-size effects. Full straight lines: crossover 
at l/^dyn ~ left ~ 0.155 (see text). 



corroborating the result (|4.25p , in spite of large correction 
terms. 

To sum up, the activation energy a(^int) is expected 
to increase monotonically with ^int all over the activated 
phase, and to interpolate smoothly between the linear 
growth (|4.23p as ^int — > ^i^t c (where the given numerical 
value of the prefactor only holds in the ^dyn — > oo limit) 
and the purely entropic limiting value ()4.25p as ^mt oo. 



4.3 Phase III (Logarithmic) 

This phase of relatively small ^int and ^dyn is illustrated 
by panel III of Figure [1] The thickness L{t) of the or- 
dered layer follows a well-defined master curve, growing 
slower than linearly with time, again with relatively small 
fluctuations between different tracks. 

Already encountered in the unhindered model [10111] , 
this phenomenon can be explained as follows. Equation 
l|4.4p shows that the application of zero-temperature dy- 
namics causes order to propagate ballistically, for ^int < 
6nt,c and <^dyn much larger than TV. When ^dyn becomes 
comparable with N, however, grains move progressively 
slowly according to their depth, with local frequencies that 
scale as (|2.3p . Writing the differential equation 



dL 
'dt 



V exp 



L 

^dyn 



we find that the thickness grows according to 



Vt 



Lit) W Cdyn In 1 + 



Cdyn 



(4.26) 



(4.27) 



and the mean jamming time for a column of size N reads 

1 I . (4.28) 



(T) « %i fexp ^ ^ 



V 



Cdyn 



This result holds all over the left of the phase diagram 
of Figure [2] (Phases I and III and the crossover between 
them). It is confirmed quantitatively by the data shown 
(empty symbols) in Figure [T] 

The laws (|4.4|) and (14. 5p are recovered for ^dyn 3> N, 
i.e., in the ballistic phase. In the logarithmic phase, when 
Cdyn N, the width of the ordered layer is predicted to 
grow logarithmically: 



L{t) « f dyn In 



Vt 

?dyn 



(4.29) 



so that the mean jamming diverges exponentially fast with 
the column size N: 



(T) 



^dyn 

V 



exp 



N 

Cdyn 



(4.30) 



4.4 Phase IV (Glassy) 

The glassy phase is found when ^int is large and ^dyn is 
small; this is by far the richest and most novel phase of 
this model. The signal for L{t), illustrated in panel IV of 
Figure [TJ is neither nearly deterministic (as in the bal- 
listic and logarithmic phases) nor totally random (as in 
the activated phase). The glassy phase corresponds to the 
'bottom' of a long column, where grain reorientations are 
at their most hindered; grains in this region are weighed 
down by those above them and, additionally, feel to the 
fullest extent the effect of the orientational frustration be- 
tween upper and lower grains. 

This phenomenon is illustrated in Figure [8] from dif- 
ferent viewpoints, using the time dependence of four ob- 
servables (for the same stochastic history which was used 
to illustrate Phase IV in Figure [1]). The jamming time 
T w 279 668 for this history is about 2.4 times larger 
than the mean jamming time for the parameters N = 70, 
6nt = 50, ^dyn = 7. The plotted observables are: 

• the thickness L{t) (see (14. 3p ) of the ordered upper layer, 

• the second component Si{t) (see p.24p ) of the pseudo- 
energy, 

• the total number i^lt) of dimers, 

• the fraction i/^ {t)/v{t) of (H — ) dimers. 

The last two quantities involve the following definitions 

of the numbers v-^ of (H — ) dimers and v |_ of ( — h) 

dimers: 



N/2 



= ^ + 0-2fe-l)(l - Cr2fc), 

fe=2 
^ N/2 

V-+ = - ^(1 - (T2k-l){l + Cr2fc), 



(4.31) 



k=2 



and of the total number v = -I- v y of dimers in a 

given configuration: 



N/2 



V- + 



-5 Ed 



<J2k-lCr2k)- 



(4.32) 



k=2 
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Fig. 8. Top to bottom: plots of the thickness L{t) of the 
ordered upper layer, the second component £i{t) of the 
pseudo-energy, the numbers v{t) of dimers, and the fraction 

i/^ (t)/v{t) of (H — ) dimers, for the history illustrating Phase 

IV in Figure [T] Full symbols: values of the observables in the 
attractor, i.e., right at the jamming time T. 



The above formulae exclude the uppermost dimer, which 
is fixed by the boundary condition (|2.12p . 

The four tracks shown in Figure [8] all show strongly 
correlated, intermittent and non- stationary fluctuations at 
all time scales, ranging from the instantaneous to scales 



of order of the jamming time (T) itself. These features 
are commonly observed in glassy systems. The existence 
of a glassy phase exhibiting this phenomenology in a one- 
dimensional model represents one of the most interesting 
outcomes of this work. 

If we examine the dynamical history depicted in Fig- 
ure in we will notice that it can be described as an alter- 
nation between two different kinds of periods: 

• Periods of quietude. Four such periods are visible in 
the figure. They are characterised by quasi-stationary sta- 
tes with a high degree of order. The thickness L{t) and 
the number vit) of dimers fluctuate around their maximal 
ground-state values of 70 and 34 respectively; the pseudo- 
energy £i{t) is correspondingly minimised, and the frac- 
tion of (H — ) dimers is either close to zero or close to unity, 
indicating a highly polarised column which is close to one 
of the crystalline attractors [/±. In some sense, it is as if 
the system has almost made its mind up to choose one of 
the two global attractors, and is dawdling in its vicinity 
with nearly no major fluctuations, during each of these 
periods of quietude. However, these long excursions do 
not in any sense anticipate the fate of the column. In the 
given example, the attractor finally chosen (full symbol) 
is close to although the system spends most of its 
time in the vicinity of the attractor [/+ with the fraction 
of (H — ) dimers typically close to unity during the periods 
of quietude. 

• Itinerant periods. During these periods of confused wan- 
dering between two consecutive periods of quietude, all 
the indicators fluctuate wildly with no particular aim in 
sight. All the observables monitored in Figure [8] are char- 
acterised by low order, with the pseudo-energy even going 
positive on occasion. 

We end up with a speculation based on a pictorial anal- 
ogy. The tracks in Figure [8] are reminiscent of those ob- 
tained in avalanche dynamics 15J, where periods of small 
random events give rise to large system-size avalanches, 
which are known to be due to stress buildup and release on 
the surface. It is interesting, using this analogy, to specu- 
late whether the itinerant periods in Figure [8] build up un- 
sustainable geometric disorder all along the column, which 
can only be relieved by a systemic choice of a nearly or- 
dered configuration (that is close to an attractor), in which 
the column then lives, until disorder strikes again in the 
form of the next itinerant period. 

Mean jamming time 

We now focus on a quantitative analysis of several as- 
pects of the glassy phase, beginning with the mean jam- 
ming time (T). Recall that in the activated phase, i.e., 
for ^int > 6nt,c and ^dyn large enough, the jamming time 
grows exponentially fast with iV (see (|4.8[) '). We now ex- 
amine the effect of decreasing ^dyn, to cross over into the 
glassy phase. 

The main effect of this is that an increasingly broad 
spectrum of local frequency scales 0„ kicks in, to slow 
down the stochastic behaviour of the activated phase. In- 
terestingly, a toy model of an assembly of independent 
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two-level systems is able to provide a clue to this crossover. 
Its details are presented in Appendix C, but the crucial 
feature is that it has two regimes - one where entropy 
dominates ('entropic'), and the other which is dominated 
by the slowest of the local frequencies ('slow'). 

The mean jamming time is in fact what would be most 
naively expected from the above competition, that is, it is 
the greater of the two times that would be generated: 

(T) ~ max (exp(a(ei„t)iV), exp(7V/edy„)) , (4.33) 

where exp(a(^int)-^) is the jamming time of (|4.8p in the 
Cdyn oo limit, while 1/0n = exp(A^/^dyn) is the slowest 
microscopic time scale of the problem. More specifically, 
this implies the following: 

• In the activated phase, when ^dyn > l/a(6nt), the result 
(|4.8|) for the mean jamming time in the ^dyn ~* oo limit is 
essentially unchanged - this corresponds to the entropic 
phase of the toy model of Appendix C. 

• In the glassy phase, i.e., for ^dyn < l/a('^int), the jam- 
ming time grows proportionally to 1/0n = exp(7V/^dyn) 
- this corresponds to the slow phase of the toy model of 
Appendix C. 

From the above, one naturally expects there to be a 
sharp transition in the fint — Cdyn plane, along the line 
defined by 

edyn(6nt) = , (4.34) 

shown qualitatively in Figure [21 For the strongest corre- 
lations, as i^int oo, the transition point ^dyn(Cint) takes 
its minimum value ^dyn(oo) = l/ooo = 2/(ln2) « 2.8854 
(see (|4.25p V On the other hand, at the boundary of the 
ballistic/logarithmic and activated regimes (^int £,^t c)' 
(|4.23p predicts a divergence of the transition point of the 
form 

^dyn(6nt) ~ 115/(6nt " 6nt,c). (4.35) 

Despite huge finite-size effects, our simulation data 
(shown as full symbols in Figure [7|) manifest the crossover 
described above. The plateau in the left part of the data 
(full horizontal line) yields the effective value Oeg ~ 0.155 
for ^int = oo and iV = 50. This effective value is very 
far from the theoretical asymptotic value (see (14.251) ). 
underlining the importance of finite-size effects. However, 
and reassuringly for our analysis, the crossover does indeed 
take place as predicted by (|4.34p . at a value of 1/^dyn ~ 
Ooff (full vertical line). 



Statistics of attractors 

We now turn to the statistics of attractors in the glassy 
phase. The question of what they are is easily addressed. 
Recall that the ground states of the model for e — 1 are 
the 2'^ configurations made up of = A^/2 — 1 dimers, 
which satisfy the boundary condition p.l2p . By construc- 
tion, these arc the possible attractors of zero-temperature 
dynamics. 



The next question, which relates to their dynamical 
attainability, is less easy to answer. A precise formula- 
tion of this question is: What is the probability Q{C) that 
the application of zero-temperature dynamics leaves the 
column in a given attractor C, starting from a uniformly 
chosen random initial configuration? Or, more physically: 
how and where does a constrained system, starting from 
random initial conditions, attain jamming? This question 
has held centre stage in theoretical [16117] and experimen- 
tal '18] explorations of granular media and many other 
complex systems, ever since Edwards postulated that the 
entropic landscape of granular systems was flat Ed- 
wards' flatness hypothesis (in the strong sense) implies 
that the attractors are sampled uniformly by the dynam- 
ics, so that Q{C) is independent of the attractor C, and 
therefore equal to the reciprocal of the total number of 
attractors. 

Our reason for introducing these issues at such a late 
stage in this paper is that the statistics of attractors are 
likely to be non-trivial only in the glassy phase. All the 
other phases indeed manifest sufficiently stochastic be- 
haviour that one would expect the entropic landscape to 
be at least approximately flat. 

A central quantity in this framework is therefore the 
dynamical entropy 

5 = -^g(C)lnQ(C). (4.36) 
c 

In the case where the attractors are sampled uniformly, 
according to Edwards' hypothesis, the dynamical entropy 
assumes its maximal value: 

5„,ax-J^ln2, (4.37) 

where = N/2 — 1. 

Measuring entropies directly via numerical simulations 
is known to be a very difficult task. Instead, we resort to 
an inspired guess. Since it seems likely that the crystalline 
attractors U± introduced in (|3.2ip will play a special role 
in the dynamics, we use them implicitly to define quanti- 
ties of interest on the attractors reached by the dynamics: 
• A global indicator is provided by the probability dis- 
tribution p{v-\ ) of the number of (-1 — ) dimers. By using 

the dimer variables rjk introduced in (I3.19P , the definitions 
(j4.3ip can be simplified as: 

^ N/2 ^ N/2 

'^+- = ^J2('^ + Vk), //-+ = -^(l-r7fe), (4.38) 

k=2 k=2 

SO that -I- V \. = V = N/2 — 1. 

If Edwards' hypothesis holds, i.e., if the 2"^ attractors 
are all equally likely to occur as attractors, the distribution 
of is binomial: 

In such a binomial distribution, the extremal values = 

and — v, corresponding to the crystalline attrac- 
tors U± (where all the dimers are of the same kind) are 
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Fig. 9. Histogram plots of the probability distribution p{h'+-) 
for A'^ = 50 (hence u = 24). Top: ^dyn ~ 10 and variable fint. 
Bottom: fint = oo and variable ^dyn- The binomial distribution 
(|4.39|) is shown as thick full lines. 



the least probable. On the other hand, if the actual dis- 
tributions obtained deviate from the binomial law (|4.39p . 
this would indicate strongly that all the attractors are not 
equally hkely, the entropic landscape is not flat, and thus 
of course that Edwards' hypothesis does not hold. 
• A local indicator of attractor structure is the correlation 
function 

{{''+- - v-+)rik) ~ 1 



Xk 



1 



(4.40) 



This correlation measures the trend for the fc-th dimer to 
be aligned with every other dimer. It vanishes identically 
if Edwards' flatness hypothesis holds. In the extreme op- 
posite situation where only the crystalline attractors U± 
are reached by the dynamics, the above correlation takes 
its maximal value Xfe = 1 for all k. 

Our numerical results for the probability distribution 

p(z/_i ) and the dimer correlation function Xk are shown 

in Figures [H] and [TUl All the data were taken for a system 
of size N ~ 50, which has v — 24: dimers that are free 
to reorient. Figure [9] shows the variation of the form of 

p{v^ ) with, first, fixed ^dyn = 10 and variable ^int, and 

next, fixed ^int = oo and variable ^dyn- Figure [TOl shows 
the variation of Xfc along the same diagnostic lines. 

In Figure [9l the binomial distribution (corresponding 
to Edwards' flatness hypothesis) is shown by a thick full 
line, with which the data for the lowest value of ^int are 
almost completely aligned. The statistics of attractors is 
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Fig. 10. Plot of the dimer correlation function Xfe against 
depth n = 2k. Parameters are as in Figure [9] 



thus very close to being uniform in the ballistic and log- 
arithmic phases, which correspond to the top layers of a 
column. As ^int increases, there is a gradual crossover to 
a non-trivial two-peaked distribution; the same trend is 
visible in the lower panel, when ^dyn decreases for infi- 
nite ^int- At the beginning of the crossover, with its small 
deviations from uniform sampling, one recognises the acti- 
vated phase, which corresponds to the middle of a column. 
By the time that the two-peaked distribution is obtained 
in both parts of Figure [9l the parameters - ^int large, and 
^dyn small - correspond clearly to the glassy phase. Here, 
attractors in the neighbourhoods of the crystalline states 
U± are eventually favoured, after long periods of systemic 
wandering. 

The above observations are reinforced by Figure [TOl 
In the upper panel, the correlation function Xk is essen- 
tially zero for low ^int, increasing progressively as ^int is 
increased. In the lower panel, i.e., for infinite ^int, the cor- 
relation function is never quite zero even for high ^dyn, 
and only manifests a stronger depth-dependence as ^dyn 
decreases. These results reinforce those found in an inde- 
pendent model which uses random graphs to model grains 
near jamming, where entropic deviations from Edwards' 
flatness occur in certain regions of parameter space |19j . 

To recapitulate, the salient feature that emerges is 
that the system prefers increasingly to live in the neigh- 
bourhood of its two global minima, the attractors C/±, as 
one goes deep into the glassy phase. As a consequence, 
the dynamical entropy decreases from a value close to its 
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maximal value (|4.37p at the boundary between Phases II 
and IV, to zero in the deepest part of the glassy phase 
(6nt ^ N, ^dyn ^ N) . Furthermore, Edwards' flatness 
is well obeyed overall in three out of four phases in this 
model; it is massively violated only deep in the glassy 
phase of this model, where the configurational landscape 
is completely rough. 

Our success in constructing this glassy phase via such 
a minimal model relied on the inclusion of two crucial 
ingredients: 

• Long-range interactions (^int > N), in order for the 
dimers to jam cooperatively, rather than independently 
of each other; 

• A broad spectrum of local frequencies (^dyn small) to 
slow down the relaxation, and thus prevent a purely acti- 
vated mechanism driven by entropy. 

It is noteworthy that both these minimal ingredients 
rely on collective effects - one to do with interactions in 
space, the other to do with degrees of freedom in time. 
This model-specific conclusion at once agrees with, and 
reinforces, general notions |13j of cooperativity in glassy 
dynamics. 

5 Discussion 

The motivation for this model came from a mental image 
of grains in a box under shaking: what could be at once 
so simple, or - as we came to see in time - so complex? 
Our first, simplest, model [8|9| involved only the effects 
of gravity on non-interacting grains - deeper grains car- 
ried the weight of grains above them, so were less free 
to move. This was modelled by using a single dynamical 
length ^dyn, representing the thickness of the dynamical 
boundary layer: grains at a depth n much less than ^dyn 
(which can be examined by setting ^dyn oo) are free to 
move, whereas those where n ~ ^dyn have lower frequen- 
cies of motion, as normal in non- Newtonian fluids. Three 
phases were found; in the 'fluidised' phase, grains flew as 
well as moved along the surface, with a relatively quick 
propagation of order down the sandbox. Grain disorder 
was essentially frozen in, in the 'glassy' phase, with a very 
slow propagation of order from the free surface. The 'in- 
termediate' phase was in some ways the most interesting, 
with a true competition between fast and slow dynamics. 

In hindsight, it is astonishing that these diverse be- 
haviours - especially the shape-dependent 'ageing' effects 
of the glassy regime ~ were manifested in a totally non- 
interacting model. A more realistic, interacting model of 
the glassy regime was presented in I10|llj . Since close- 
packed grains can typically not diffuse spatially, it was 
sufficient to model a column, rather than a box, of grains. 
In the model of ^10111] , grain motions were constrained not 
just by the masses, but also by the orientations of grains 
above them, thus generating directional long-range inter- 
actions. The effect of compaction around jamming was 
modelled by a single local field h„, representing the ex- 
cess void space [20j for grain n, which could be minimised 
by a suitable choice of grain orientation. Additionally, in 
this model, grains were allowed to have arbitrary shapes. 



so that the the disordered orientation of a grain could 
occupy any volume e, and, correspondingly generate any 
void space. The propagation of order in this model pro- 
ceeded from the free surface to the base, and was 'causal in 
space' - in that while upper grains constraint lower ones, 
the converse was not true. 

It took us some time and several explorations to realise 
that while the model of |10|llj had at least the flavour of 
the interactions needed to model a jammed glassy phase - 
e.g. the constraining effect of long-range grain correlations 

- its lack of slow dynamics (except those arising from the 
trivial effect of grain masses) was a direct result of its spa- 
tial causality. Essentially, provided a grain was not blocked 
down by the weight of other grains, it was free to orient 
itself subject only to the orientations of grains above it- 
self - that is, we were modelling the behaviour of the top 
layers of a jammed column of grains, which never felt the 
undertow of the base. It was small surprise, therefore, that 
the ordering dynamics for ^dyn oo were ballistic. 

To model a column of grains with spatial inhomogenei- 
ties - that is, a column with a top, a middle, and a bottom 

- we discovered that orientational constraints needed to be 
inserted in a non-directed way. This enabled us to model 
frustration - in this context, the need of a given grain to 
balance the effects of two competing local fields /i„ and jn 

- which led in its turn, to slow dynamics. Still keeping the 
orientational constraints of previous models [8 9 1 0111] via 
the field /i„, as well as the effect of 'gravitational slowing 
down' via ^dym we therefore introduced here the notion 
that grains were also constrained by grains below them, 
via the field j„ which propagated over a correlation length 

Cint- 

From this very heuristic and pictorial modelling has 
emerged a model column of grains that manifests all the 
complexity of earlier models |10|llj . and adds some more 
via the introduction of the activated and glassy phases. 
Our main success is of course in the realisation of a glassy 
phase for which a minimal combination of two physical 
ingredients - strong, bi-directional, orientational correla- 
tions and a broad spectrum of local frequencies arising 
from a natural depth- dependence - appears to be neces- 
sary. The richness of this phase is worthy of further explo- 
ration, especially to do with issues concerning higher-order 
correlations and ageing. 

Finally, we mention that our model provides some ra- 
ther interesting and general insights into the nature of op- 
timisation - the granular column modelled here reaches its 
ground states in strikingly different ways, in the four dy- 
namical phases mentioned above. It is tempting to think 
of these phases as representing different spatial parts - 
'top', 'middle' and 'bottom' - of a column, and to con- 
nect their different routes to compaction with the issue 
of inhomogeneities in real granular media [1]. While this 
picture is an appealing one, one should remember that 
the four phases of this model were obtained by varying 
^int and fdyn; the translation of our results to apply to a 
real column would involve the natural apparition of such 
variations as a function of depth. 
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Appendix A. Exact results for = 4 and = 6 

In this Appendix we derive analytical results on the zero- 
temperature dynamics of the model investigated in the 
body of the paper, for small systems. We concentrate 
onto the statistics of attractors, i.e., the probability that 
the system is absorbed in each attractor, starting from a 
random initial configuration. Exact results will be succes- 
sively obtained for systems of sizes = 4 and N = 6, for 
arbitrary ^int and ^dyn- 



The case N = A 

The column made of 4 grains with the boundary condition 
(|2.12p has three free orientations : (T2 ± 1, 0-3 = ±1, 0-4 ^ 
±1, and therefore eight configurations, in lexicographical 
order: 



Ci = 



+ +: C2 = 



+ +-, C3 = + 

- ++, Ce = + 



-+-, (A.l) 



The two attractors are the dimerised configurations 
Cq and C7, respectively corresponding to 772 = +1 and 
?72 = —1. Our main goal is to determine the probabilities 
Q^^^ and Q^'^^ that the system is absorbed in each of these 
attractors, starting from a random initial configuration. 
The relevant information is encoded in the non-uniformity 
parameter 

^ = (??2> = Q('^ - (A.2) 

The zero-temperature dynamics consists of a certain 
number of moves between configurations. For instance, Ci 
may be updated into the following configurations at the 
following rates: 



/I 03 /I 04, 

Ci > C3, Li > C2. 



(A.3) 



As announced in Section 3.3, this dynamics is independent 
of ^int for A^ = 4. It can be represented as an 8 x 8 Markov 
matrix M, such that the occupation probabilities Pa{t) 
of the configurations Ca (a = 1, . . . , 8) obey the forward 
Kolmogorov equation [211221 



(A.4) 



The absorption probabilities of the attractors can be 
derived by means of the following approach. Let Qa be 
the probability of being eventually absorbed by configu- 
ration Cc, starting from the initial configuration Ca- For 

(c) 

a fixed configuration Cc, the absorption probabilities Qa 
obey the backward Kolmogorov equation [21l22j 



(A.5) 



1. For a 



complemented by the boundary condition Q 
random initial configuration, we have therefore 



By solving (|A.5[) successively for both attractors [c — Q 
and c = 7), we are left with the following explicit result, 
for arbitrary rates 0„: 



0304(03 - 04) 



2(02 + 03)(02 



l)(03 



■04 



(A.7) 



In the present model, where 0„ = x^^^ (see (I2.3p ). this 
result simplifies to 



A = 



(1 a;dyn)2:dyn 

2(l + Xdy„)'(l + 2;2yJ- 



(A., 



The non-uniformity parameter Z\ is a small positive quan- 
tity, meaning that the crystalline attractor Cg = is al- 
ways slightly favoured by the dynamics. It coincides with 
Ai plotted in the upper panel of Figure [TTl It vanishes in 
both limits Xdyn (i.e., ^dyn 0) and 1 (i.e., 

Cdyn — > oo). Its maximum A w 0.012465 is reached for 
Xdyn ~ 0.6253. 

The case A^ = 6 

The case of a column made of A^ = 6 grains can still be 
dealt with by analytical means, although the final expres- 
sions are much lengthier. The system has 32 configura- 
tions. Its four attractors are the following configurations 
(relabelled as a = 1, . . . , 4 for convenience): 



C3 



C2 
C4 



(A.9) 



The relevant information is encoded in the three non- 
uniformity parameters 

^i = (r/2) = Q(i) + Q(')-g(^)-g(^\ 

A2 = (%) = Q(l) - Q(2) + g(3) _ g(4)^ (A.IQ) 

z\3 = (772??3) = g(^)-g(^)-g(^) + g(^). 

As announced in Section 3.3, A^ = 6 is the smallest 
system size such that the dynamics depends in a non- 
trivial way on the parameter ^int- The Markov matrix M 
assumes two different expressions for < Xjnt < 4> and 
< a^int < 1, where the inverse golden mean is given in 
(I3.15p . For each of these two phases, the 32 x 32 Markov 
matrix M has been generated, and the backward equa- 
tions (IA.5P corresponding to each of the four attractors 
have been solved analytically^ We thus obtain the follow- 
ing expressions for the non- uniformity parameters, in the 
range < x-mt < (t>- 



Ai^ 



(1 Xdyn)a;dyn 

2(l + rrdy„)'(H-a:^y„)' 



^2 = {X%^ - l)(a^dyn + Xlyn + 1) 

^ (^i ,^3 ,„2 ,1^™6 -P24(xdyn) fAll) 
^ {^dyn + ^dyn + ^dyn + '^l^dyn j^^^^ ^' V^'-^-^j 



As^ -(Xdyn-l?{xly,, 



i)(^: 



X {x'^ 



dyn 



''dyn 



dyn ^>dyn ^(^^^^) ' 



dyn ~ dyn 
5 f25(a;dyn) 



1) 



with help of the software MACSYMA. 
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and in the range (j) < Xint < 1: 



^1 = - 



^2 - {4,n 



'^^'^ D{xdyn) ' 
l)(4yn + 



'^dyn 



^ \-^dyn ' -^dyn ~ "^dyn ^ -^y-^dyn ]J(^x^ ) ' 



(A.12) 



,3 

dyn 



^3 = (l-4yn)(-T^ 

^ (^^dyn + 2^dyn + 



^dyn 



''dyn 



s ^30(2;dyn) 



£'(a;dyn) 

where we have introduced the following polynomials: 



J2 



X (a;** - 



(a;^ + a; + 



- l)2(a;3 



l)(a;2 



- X 
^3 



l)(a;4 + 



- X 



+ n 



a;2 + l), 
(A.13) 

434x1'^ 
- 533x9 
x^ 



+ 5x23 + 10x^2 _^ 152;21 _|_ 92,20 

51x1^ - 104x1^ - 203x^6 - 308x^5 - 
515x13 - 605x^2 - 634x" - 609x1" , 
449x* - 330x^ - 218x6 - ISOx^ - 7] 
27x3- 3x2 +X + 1, 

(A.14) 

,25 

+ 335x^0 + 547x19 + 842xi** + 1177x" + 1556x1^ 
+ 1900x15 + 2193x1^^ + 2349x13 + 2356x1^ 
+ 2223x11 + 1955x1" + 1609x9 + 1242x« + 888x^ 
+ 590x6 3572,5 207x'^ + 101x3 + 43x2 
+ 17x + 4, 



P25(x) = 6x25 + 15x24 + 40x23 + 95x22 



- 292x 
2283x22 



,26 



= 6x30 _^ 21x29 + 54x28 + 140x27 
+ 546x25 + 949x24 + 1531x23 

+ 6457x18 
+ 8236x14 

+ 7793x13 + 7011x12 + 5978x11 _,_ 4852a;io 
+ 3730x9 _^ 2682x8 + 1811x7 + 1142x6 
+ 343x4 + 160x3 + 64x2 + 20x + 4, 



(A.15) 



21 _ 


h 4277x2" - 


1- 5377x19 


17 _ 


h 8025x16 - 


h 8316x15 


13 _ 


1-7011x12- 


h 5978x11 



p- 



38 (a;) = 4x38 + 28x37 + 109x36 + 341x35 _^ 9092;34 

+ 2114x33 + 44 1 6x32 + 84 24x31 i4836a;30 



662x5 
(A.16) 



'37105 



„28 



53309x27 ^ 72124X 



,26 



+ 24262x29 

+ 92 1 05x25 + 111177x24 + 12 7 001x23 
-t- 137 1 90x22 + 140 74x21 -|- 135038x2" 
+ 122612x19 + 104406x18 + 82925x17 - 
+ 40790x15 + 24131x14 + 11829xi3 -|- J 
- 720x11 - 2614x1° - 2891x9 - 2404x8 - 1693x7 



60963x16 

■nn™12 



1031x6 - 555x5 
llx - 2. 



269a 



112x3 _ 39^,2 



(A.17) 

The non-uniformity parameters Ai are plotted in Fig- 
ure [11] (note the powers of 10 in the vertical scales). For 
< Xint < (t> (upper panel), the Ai are typically small. 
They vanish in both limits Xdyn — *■ and Xdyn — *■ 1- Ai is 
always positive, and it coincides with A of the case = 4 
(see (jA.SP ): A2 turns from negative to positive for Xdyn ~ 
0.2598, whereas A3 is always negative. For cj) < Xint < 1 
(lower panel), the Ai are typically larger. Only A2 is the 
same in both phases. The Ai again vanish in both lim- 
its Xdyn and Xdyn 1, except for Ai — —5/96 for 
a^dyn = 1 and ^3 = 1/2 for Xdyn — 0. Finally, Ai turns 




dyn 



o A xlO 
• aIxIOOO 




Fig. 11. Plots of the non-uniformity parameters Ai (i = 
1,2,3) characterising the attractor statistics of a column of 
iV = 6 grains. Top: < xint < (see (|Xll|). Bottom: 
(j) < a;int < 1 (see (|A.12[) \ Note the powers of 10 in the vertical 
scales. 



from positive to negative for Xdyn 
is always positive. 



0.6581, whereas A3 



Appendix B. Biased Brownian motion on an 
interval 

In this appendix we consider biased Brownian motion on 
the interval < x < L, characterised by its drift velocity 
V and its diffusion coefficient D. The endpoint at x = 
is reflecting, whereas the endpoint at x = L is absorbing. 

This mixed type of boundary conditions is referred to 
as the transmission mode [21j . It ensures that, starting 
at a given position x = £, the particle hits the absorb- 
ing endpoint with certainty in some finite random time 
T, referred to as the absorption time. Our main goal is to 
determine the distribution p(T) of this absorption time, 
and especially its first few moments. We shall mostly con- 
centrate onto the £ — > limit, where the particle starts 
from the immediate vicinity of the reflecting endpoint. 

Let P{x,t) be the probability density for the parti- 
cle to be at point x at time t, and J(x, t) the associated 
probability current. We have 



dP 
'dt 



dJ_ 

dx 



0, J = VP-D 



dP 

dx ' 



(B.l) 
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with the initial and boundary conditions 

P{x,0)^5{x-i), J(0,t) = = 0. (B.2) 
Introducing the Laplace transform 

/■OO 

Pix,s) = / P{x,t)e-''* dt, (B.3) 

"'0 



the above equations imply 

dP „ d^P 



sP + V- 



D- 



dx dx^ 
with boundary conditions 



5{x - £), (B.4) 



VP{Q,s)^D^^^^^^P{L,s) = Q. (B.5) 
Equation (jB.4p implies the matching conditions 

p{i+,s)^p{i-,s), 

dP{£+,s) _ dP{£-,s) _ (B.6) 
dx dx D 

Consider first the homogeneous equation obtained by 
setting the right-hand side of (jB.4[) equal to zero. Looking 
for a solution to this equation at fixed s of the form P = 
e''^, we obtain the quadratic equation s + Vr — Dr^ = 0, 
whose two roots read 

V-W V + W , n ,1/9 



(B.7) 



2D ' 2D 
The solution obeying the boundary conditions (jB.5[) reads 

A{ri e"^'-" - r2 e''^^) (0 < a; < £), 



P{x,s) 



(B. 



Finally, the constants A and B are determined from the 
matching conditions (jB.6p : 



A 



B 



r2 e — ri e 



D{r2-n){r2 e-'-i^-ne-'-^^)' 
The survival probability of the particle at time t, 

S{t) = [ P{x,t)dx, (B.IO) 

is nothing but the probability that the absorption time T 
is larger than t: 

/oo 
p{T)dT. (B.ll) 

We have therefore, in Laplace space 



S{s) = / P{x, s) dx, p{s) = (e'"^) = 1 - sS{s). 
Jo 

(B.12) 



The solution (|B.8p . (|B.9p yields after some algebra 



and 



r2(e-'-i^ 



-riL 



) -ri(e- 



r2 e-'''^ - n e"''^^ ' 



(B.13) 



(B.14) 



In the following we restrict the analysis to the limiting 



£-> 0, 



(B.15) 



where the particle starts from the immediate vicinity of 
the reflecting endpoint. The result (|B.14p simplifies to 



Pis) 



r2 - n 



r2 e^''^^ - ri e""^^-^ ' 



(B.16) 



i.e., explicitly, 

m 



2We 



VL/{2D) 



(^W + V) e^^/(2^) + {W-V) e-'^^/f^^) ' 

(B.17) 

where W has been defined in (|B.7p . 

The moments of the absorption time T can be derived 
by expanding the result (jB.17p as a power series in s, as 
p(s) = 1 - {T)s + (T2)s2/2 + • • • We thus obtain 



(T) = -^iVL -D + Dc-^^/^), 



V 

(T2) = -^{V^L'^ - 41)2 ^ 2D{WL + D)c-^^/^ 
+ 2D\-^^^'D), 

(B.18) 

so that the reduced variance of the absorption time T, 

(T)2 (T)2 ^ ' 



reads 



Kt = D 



2VL -5D + A{VL + D)e-^^'^ + Dc^^vl/zj 



(V^i _ D + 75e-^^/^)2 

(B.20) 

The above results have different kinds of behaviour in the 
following cases. 

• Ballistic phase {V > 0). In this case, the drift brings the 
particle toward the absorbing endpoint. The mean absorp- 
tion time 

(^) - f - ^> (B.21) 

grows linearly with L, according to a ballistic law with 
velocity V . Note the negative correction due to diffusion. 
Fluctuations of the absorption time around its mean are 
asymptotically Gaussian. Their reduced variance, 



2D 

vT' 



(B.22) 



falls off as 1/L. 



20 



J.M. Luck, Anita Mehta: Dynamical diversity and metastability in a hindered granular column near jamming 



• Diffusive point {V = 0). This special case corresponds 
to pure diffusion. The expression (jB.17|) simplifies to 

1 



(B.23) 



(B.24) 



coah{{s /D)^/^Ly 
The mean absorption time 

(T) - To = 

^ ' ^ 2D 

defines the diffusive time scale Tq, which grows quadrati- 
cally with L. The Laplace transform can be inverted ex- 
plicitly in this case. The dimensionless ratio r = T/Tq has 
a non-trivial distribution: 

oo 

Pir) = f 5](-l)'=(2fc+l)exp(-(2fc+l)27rV8), (B.25) 



k=0 



with moments (r) = 1 (by construction), (r^ 
(t^) = 61/15, and so on. We have thus 



1 = 1 
3 



5/3, 



(B.26) 



• Activated phase (V < 0). In this case, the drift brings 
the particle toward the reflecting endpoint. It is therefore 
very improbable that the particle sits by chance near the 
absorbing endpoint. The absorption mechanism is there- 
fore activated. The mean absorption time, 



(T) 



ZJel^l^/^ 

y2 ' 



(B.27) 



is found to grow exponentially with the length L. The 
corresponding activation energy per unit length reads 



D ' 



(B.28) 



The distribution of the absorption time can be checked to 
be asymptotically exponential. In particular, the reduced 
variance, 



1 



2iVL + 3D) 
D ' 



\V\L/D 



(B.29) 



converges exponentially fast to the limiting value unity, 
characteristic of the exponential distribution. 
• Critical phase {V small, L large). For a large length L, 
the distribution of the absorption time exhibits a finite- 
size scaling form in a narrow interval of V around the 
diffusive point V = 0, whose width scales as 1/L, where 
the dynamics interpolates between the ballistic and the 
activated phases. Let us introduce the dimensionless finite- 
size scaling variable 



VL 



(B.30) 



which is twice the Peclet number introduced e.g. in ^21^ 
The moments of the absorption time scale as 



(T> = To 
(T2) = T 



2(z- l + e"^) 

2 4(z2 + 8(3z -I- 1) e^^ -I- 8 e'^^ 



(B.31) 



The reduced variance therefore depends continuously on 
z according to 



2z-5 + 4(z + l)e-^ 



(B.32) 



Appendix C. Mean hitting time for two-level 
systems 

In this appendix we consider an assembly of N indepen- 
dent, albeit not identical two-level systems, described as 
spins s„ = ±1 (n = 1, . . . ,N). The spins are flipped ac- 
cording to independent Markov processes whose rates 



w{Sn = +1 ^ S„ = -1) = w{Sn = -1 



= +!) = 



(c.i) 

depend on the label n in an arbitrary fashion. The sta- 
tionary state of this dynamical process is an equilibrium 
state where the 2^ configurations C = {si, . . . , sat} are 
equally probable. The above dynamics indeed obeys de- 
tailed balance with respect to the uniform measure. 

Let Ct denote the configuration of the system at time 
i, starting from a given random initial configuration Co, 
For a given stochastic history of the system, we introduce 
the hitting time 



T = min{i | Ct = C*}, 



(C.2) 



defined as the first time the system visits a given configu- 
ration of reference, say = {+1, . . . , +!}• 

It will be sufficient for our purpose to evaluate the 
mean hitting time (T), averaged both over the random 
initial configuration Cq, chosen with the uniform equilib- 
rium measure, and over the stochastic dynamics. Along 
the lines of the introduction of [21], it is advantageous to 
consider simultaneously the probability 



Pc.Ca{t)^Proh{Ct^C\Co} 



(C.3) 



for the system to be in configuration C at time t, knowing 
that it was in configuration Co at time 0, and the proba- 
bility 

Fc Co (t) = Prob{Ct =C, C't^C for all < i' < t | Co} 

(C.4) 

for the system to be in configuration C for the first time 
at time t, again knowing that it was in configuration Co 
at time 0. These two quantities are related by the integral 
equation 



Pc.Coit) 



Fcfio{t')Pc,c{t~t')dt'- 



(C.5) 



Let us average this formula with respect to the initial con- 
figuration Co, chosen with the uniform equilibrium mea- 
sure. The left-hand side equals 1/2^, for all configurations 
C and all times t > 0. The return probability Pc,c{t — t') = 
R{t~t') is also independent of C. As a consequence, the av- 
erage over Co of Fcfia (^') defines some average first-passage 
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probability F{t'), which does not depend on C either. The 
return probabihty R{t) and the first-passage probabihty 
F{t) obey the convolution equation 



f Fit')R{t-t')dt', 
^ Jo 



i.e., in Laplace space, with the notation (|B.3|) . 



(C.6) 



(C.7) 



In the present case of iV independent spins, the return 
probability factorises as 



TV 



(C.8) 



where 



r„(t) = Prob{s„(t) = s„(0)} = -(1 + (s„(t)s„(0))). 

(C.9) 

The temporal correlation function of each spin variable 
exhibits a pure exponential decay at equilibrium: 



(s„(i)s„(0)) =e 



SO that 



N 



(C.IO) 



(C.ll) 



The return probability tends toward the limit R{oo) = 
1/2^, as it should. Its Laplace transform therefore has 
the following behaviour as s — > 0: 



(C.12) 



where the finite part C is given by the convergent integral 

C = j f [](l + c-2«"*)-l j di. (C.13) 

The mean hitting time under consideration is just the 
mean value of the first-passage time: 



f°° dF 
(T)= tF{t)dt = - — {s = Q). 



(C.14) 



Expanding \C7\ to first order in s, we obtain (T) = C, 
hence our final result: 



(T) 



X (n(l + e-^^"*)-lj dt. (C.15) 



Expanding the product and integrating the exponentials 
term by term, we obtain 



(C.16) 



In this expression / runs over the 2^—1 non-empty subsets 
of {1, ... , N}. We thus have for the first few values of N: 



^ = 2: (T) = U- + - 



N ^3 



1 



2X01 02 01 -1- 02/' 

1/1 1 1 

T = - — + — 

2V01 02 03 

1 1 

H \ 

01 +02 01 + 03 

1 1 

+ 



(C.17) 



02 + 03 01 + 02 + 03 

The situation of interest in the body of this paper is 
where the fiipping rates read 0„ = x^y^ (see (12.31) ). where 
.Tdyn assumes any value in the range < x^yn < 1- The 
following special cases can be worked out more explicitly. 
• Uniform case {x^yn = !)• In this case, all the flipping 
rates are equal (0„ = 1). Equations (|C.15|) and (|C.16|) 
simplify to 



(r) = /°°((i + e-^r-i) dt^l 



(C.18) 

where m is nothing but the number of elements of the set 
/ entering (|C.16p . We thus obtain the estimate 



(T) « — . 



(C.19) 



• Binary case (xdyn = 1/2). The problem also simplifies in 
this case, where 0„ = 1/2". Consider indeed (|C.16p . The 
denominator of the generic term can be recast as 



(C.20) 



nel 



nel 



The sum in the right-hand side is nothing but the binary- 
digit expansion of a generic integer M = 1, . . . , 2^ — 1. 
Therefore 

2"-l 



We thus obtain the estimate 



(T) 



M=l 



A^ln2 



(C.22) 



This particular value demarcates two different kinds of 
behaviour. The change of variable C = ln(2t)/A in (|C.15|) . 
with the notation A = jlnxdynl, indeed yields 



A 



N 



= i y ^ n (l + exp(-e(«-")^)) - 1 j eC^ dC. 

(C.23) 

The factors in the product are approximately equal to 2 
when the difference C. — n is very negative, and to 1 when 
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it is very positive. As 
accuracy, we obtain 



a consequence, with exponential 



max (2 



N 



dyn 



dyn) 

(1/2 < Xdyn < 1), 
(0 < Xdyn < 1/2). 



(C.24) 



These estimates are exact up to Xdyn-dependent prefac- 
tors, as we now successively show for 1/2 < Xdyn < 1 and 
for < Xdyn < 1/2. 

• Entropic phase (1/2 < Xdyn < 1)- In this phase, where 
the rates 0„ exhibit a rather mild dependence on n, the 
exponential estimate (T) ~ 2^ has an entropic origin. 
Just as in ()4.24p . it scales as the ratio between the initial 
phase-space volume (all the 2^ configurations) and the 
final one (one single configuration of reference) . 

A more accurate estimate of (T) goes as follows. Re- 
writing each factor of the product entering (jC.lSp as 



1 + e ^*^dy„ ^ 2 e~*^dyn cosh(tedyn) 



we obtain 
with 



(T) « 2^M(a;dyn), 



(C.25) 
(C.26) 



-tXdya/ (l-aJdyn) 



Y[ cosh(tedy„) dt, 



n=l 



(C.27) 

where both the product and the integral are convergent. 
The series expansion of the amplitude yl(a;dyn) as Xdyn 
1 can be derived by expanding the infinite product as a 
power series in t and integrating term by term. We thus 
obtain 



^(a;dyn) = (1 - 2^dyn) + -^C^ - Xdyn)'^ 

5 19 

+ 2(l-2:dyn)' + ^(l-a:dyn)^ + 



(C.28) 



The amplitude A{xdyn) vanishes linearly as a;dyn — * 1, so 
that ([Cl^ and ((C?^ are compatible. The result (fC^ 
suggests that ^(a;dyn) diverges linearly as Xdyn ^ 1/2+: 



^(s^dyn) 



In 2 



4(a;dy„ - 1/2) ' 



(C.29) 



Finally, we have A{xdyn) = 1 for Xdyn ~ 0.6396. 
• Slow phase (0 < Xdyn < 1/2). In this phase, the scale of 
the hitting time is given by the slowest time scale of the 
system: (T) ^ l/0jv = 1/x^^^. 

A more accurate estimate of (T) can be derived as 
follows. Setting n — N — j in (|C.16|) . we obtain 



(T) 



^(a;dyn) 



''dyn 



where 



^(a^dyn) = IY1 




Fig. 12. Plot of (l/iV)ln(r) against |lna;dyn|, illustrating the 
behaviour of (T) in the various phases. Full lines: Data for 
N = 6, 10, 14, 18, and 22. Dashed straight lines: exponential 
estimates (IC.24[I . Full symbols: values of Xdyn where A{xdyn) = 
1 or B{xdyn) = 1 (see text). 



and where J runs over all the non-empty subsets of the 
integers {0, 1,2,.. .}, i.e., explicitly 



^(a^dyn) = + 2Xdyn + ^X^^n + ^^ly^ + H 

+ 254y„ + 44x6y„ + 94x^y„ + . 



''dyn 



(C.32) 

The result (IC.22[) suggests that the amplitude i3(a;dyn) 
diverges linearly as Xdyn l/2~: 



B{Xdyn) 



In 2 



4(1/2 - Xdyn) 



(C.33) 



Finally, we have -B(xdyn) = 1 for 2^dyn ~ 0.2662. 

Figure [12] shows a plot of (1/A^) ln(r) against |lnxdyn|, 
for system sizes N ranging from 6 to 22. The data exhibit 
a sharper and sharper crossover between both exponential 
estimates (|C.24|) . shown as dashed lines. The data for all 
the values of N intersect the dashed lines very near the 
theoretical values Xdyn ~ 0.2662 and Xdyn ~ 0.6396, where 
the amplitudes i?(xdyn) and A{xdyn) are equal to unity. 
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